subroutine vmd_ic(f0vec, dfu, uxMat, E0C, phi0C, phi0T, fvc0T, t)    ! outputs

  use dimensions

  use fft_wrap, only   : fft_forward, fft_backward

  use parameters, only : phi0, kz, vPhase, u, &                       ! real,     phys
                         ny, iWaveFrame, iModx, iModz, &              ! integer,  phys
                         ModAmp,  &                                   ! real,     phys
                         restart, &                                   ! integer,  num
                         idfMod, &                                    ! integer,  num 
                         pi, iunit, isign                             ! constants


  implicit none

!-- interface input/output --

  real, intent(out)    :: f0vec(jv)
  real, intent(out)    :: dfu(jz, jv, jx, jDirection)
  real, intent(out)    :: uxMat(jx, jz, jDirection)
  complex, intent(out) :: E0C(jz)
  complex, intent(out) :: Phi0C(jz)                   ! restart value overwriten
  complex, intent(out) :: phi0T
  complex, intent(out) :: fvc0T(jv)
  real, intent(out)    :: t



!-- local variables ---

  integer              :: ix, iz, iv, iDirection  ! counters in r-space  

  real                 :: MoDinZ(jz), MoDinX(jx) 
  integer              :: jzOLD, jvOLD, nyOLD

  real                 :: v, f0
  real                 :: meandf(jv)

  complex              :: fvc0(jv)



  f0(v)=( 1/sqrt(2*pi) )*exp(-.5*(v+vPhase*iWaveFrame)*(v+vPhase*iWaveFrame))

!----------------------------------------------------------------------------------
!-- create vector "f0vec" ---------------------------------------------------------


      do iv=1,jv
        v = vvec(iv)
        f0vec(iv) = f0(v)   
      end do


      open(1,file='stuff_vs_v_initial.txt')
      write(1,*)'% v  f0  <f>'
      meandf=0
      do iv=1,jv
         do iDirection = 1,2
            do iz=1,jz
               do ix=1,jx
                  meandf(iv) = meandf(iv) + dfu(iz,iv,ix,iDirection)
               enddo
            enddo
         enddo
         meandf(iv)=meandf(iv)/(jx*jz)
         write(1,501) vvec(iv), f0vec(iv), f0vec(iv) + meandf(iv)
      end do
      close (1)


      open(13,file='f0.txt')
      write(13,*) '% V/Ve  f0Maxwell  f0Mod'; 
      do iv=1,jv
         write(13,501) vvec(iv), f0(vvec(iv)), f0vec(iv)
      enddo

501   format(  6400(e12.6,2x)  )



!----------------------------------------------------------------------------------
  open(5, file='output', position='append')

  if(restart.ne.0) then         ! read restart info

      open(4,file="restart")

      write(*,*)
      write(*,*)'reading from restart file'
      read(4,*) jzOLD, jvOLD, t, nyOLD
      if (ny.ne.nyOLD) then
         write(*,*)'from restart file, ny.ne.nyOLD, kill program'; stop
      endif
      if(jvOLD.ne.jv)then
         write(*,*)"from restart file, jvOLD.ne.jv, stop execution and return"; stop
      end if
      if(jzOLD.ne.jz)then
         write(*,*)"from restart file, jzOLD.ne.jz, stop execution and return"; stop
      end if
          
      read(4,*) E0C, phi0C, phi0T  
      read(4,*) uxMat  
      read(4,*) dfu  

      close(4)      

      write(*,*)'FINISHED reading from restart file' ; write(*,*)                  

   else                   ! start fresh
       t = 0
       dfu=0.0
       phi0T = phi0
       do iz=1,jz
          E0C(iz) =  -iunit * kz * phi0 *exp ( iunit * kz * (iz-1) * dz  )
       end do
       write(*,*)' fresh start, jz=',jz,'   jv=',jv,'    ny=',ny
 
       do iDirection = 1, 2
          uxMat(:,:,iDirection) = u * iSign(iDirection)
       enddo
      
       if(idfMod==1)then

          if( iModx.ne.1 ) iModx=0 ;   if( iModz.ne.1) iModz=0
         
          MoDinZ(:) =    iModz * cos(  kz * zvec(:) )
          MoDinX(:) =    iModx * cos( dkx * xvec(:) )
         

          do iDirection = 1,2
             do iv=1,jv                                                      
                do ix=1,jx
                   do iz=1,jz
                      dfu(iz,iv,ix,iDirection) =  ModAmp * f0vec(iv) *  (  MoDinZ(iz) + MoDinX(ix) )
                   enddo
                enddo
             enddo
          enddo !   iv and iDirection
                 
          write(5,*)

          if( (iModx==1).or.(iModz==1) ) then
             if(iModx==1)   write(5,*)' df is  modulated at wavenumber kx, kx/kz =',dkx/kz,'   with amplitude =', ModAmp
             write(5,*) 
             if(iModz==1)   write(5,*)' df is  modulated at wavenumber kz','   with amplitude =', ModAmp
             write(5,*)
          endif
       endif !  idfMod==1
    end if  !  end restart test

!----------------------------------------------------------------------------------

close(5)


      !-- creating complex vector out of real and its FFT. Why? --- 
                  
      fvc0(:)=f0vec(:)      
      call fft_forward(jv,fvc0,fvc0T)   !  fvc=FT(dfu)
      fvc0T=fvc0T/float(jv)             !  normalized velocity  FT of equilibriium
      
 
      call fft_forward(jz,E0C,E0C)
      E0C=E0C/float(jz)
      call dealias(jz,E0C)   ! in Fourier space

      Phi0C(:) = iunit * E0C(:) / kZvec(:)

      call fft_backward(jz,E0C,E0C)

      call fft_backward(jz,Phi0C,Phi0C)  ! back in x-space









end subroutine vmd_ic
